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Abstract. I consider the problem of deriving couplings of a statistical model from mea- 
sured correlations, a task which generalizes the well-known inverse Ising problem. After 
reminding that such problem can be mapped on the one of expressing the entropy of a 
system as a function of its corresponding observables, I show the conditions under which 
this can be done without resorting to iterative algorithms. I find that inverse problems 
are local (the inverse Fisher information is sparse) whenever the corresponding models 
have a factorized form, and the entropy can be split in a sum of small cluster contribu- 
tions. I illustrate these ideas through two examples (the Ising model on a tree and the 
one-dimensional periodic chain with arbitrary order interaction) and support the results 
with numerical simulations. The extension of these methods to more general scenarios is 
finally discussed. 



1. Introduction 

The estimation of the direct interactions among the microscopic constituents of an ex- 
tended system is a problem that has recently received considerable interest from the liter- 
ature of several communities (biology [HE], genetics [31 E], neuroscience [3 El [7], economy 
El [ID]). This is mainly due to the availability of large datasets across various fields, 
which have introduced the possibility to directly fit from data the interaction structure 
(e.g., the wiring pattern of neurons, the regulatory network of genes, the cross-influence 
of traders) of a system which is operating according to an unknown, or a partially known, 
mechanism. Despite this abundance of data, the problem of estimating the most appro- 
priate model in order to describe a system is far from being solved: even if the maximum 
entropy principle (MEP) is invoked in order deal with the issue of model selection 112]. 
one is typically left with a difficult inference problem. In particular fitting a model un- 
der the MEP requires matching a set of observables estimated from data with the ones 
predicted by the model (the so-called inverse problem), a task which is known to be gener- 
ically hard [13]. While the Boltzmann learning algorithm can be used to obtain accurate 
solutions for the inverse problem only by investing a large amount of computational power 
[H], mean-field techniques [US EH EEZ1 E3 EH EI] , message-passing methods [2T| 1221 1231 123] 
and pseudo- likelihood approximations |25] [27>] have proved to be very effective in order to 
fit approximately the interactions of the system in a short amount of time. 
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In this work I show that for a specific class of models the inverse problem can be solved 
exactly and efficiently. In particular, if the probability distribution of the system has a 
factorized form, the entropy of the system can be written as a sum of cluster contributions 
which are easy to manipulate analytically. This is consistent with the results of \27\ [28] , 
who show that efficient, fast inference procedures can be constructed if the entropy of the 
system can be approximated by the sum of a suitable set of small cluster contributions. 
For these systems one can see that the robustness of the inverse problem is linked with 
the factorization property of the probability distribution (equivalently, the additivity of 
the entropy). This is a consequence of a general fact: inverse problems are typically more 
stable with respect to the direct ones because the response of couplings due to a change in 
the observables is much more localized than the response of the observables due to a shift 
of the couplings. 

In the following section I introduce a very generic setup for the inverse problem, in or- 
der to stress the generality of the results. In section [3] its exact solution is presented for 
two specific cases. In [4] the results are scrutinized against numerical simulations involving 
synthetic data, and in [5] the method introduced is extended to more general cases. In 
section [6] some conclusions are drawn. 



2. The inverse problem 

I consider the problem of estimating a set of M couplings g = (fli, ... ,5m) describing a 
probability density for a set of N binary spins s = (s±, . . . , sn) G { — 1, 1}^ of the form 



where <j> = (<pi, . . . , 4>m) is a se t of known functions of s (usually referred as potentials, or 
sufficient statistics), and the partition function Z{g) enforces the normalization of p(s). 
In the following discussion we will assume all of the 0«(s) to be monomials, i.e., to be 
functions of the form 



where is a given subset (cluster) of the vertex set V = {1, . . . N} associated with the 
operator fj,. Within this general formulation, the Ising (or graphical) model can for exam- 
ple be constructed by using single-body clusters {i}^ and two-body clusters {i,j}f<j=ii 
together with their associated couplings usually denoted with {hi}f =l and { Jij}^. J=1 . Fig- 
ure [T] shows factor graph associated with the probability density ([TJ in a generic case. In 
order to estimate such model, a set of T i.i.d. observations s = {s^^}f =1 oi the system is 



(1) 




(2) 



= n 
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Figure 1 . Example of a typical factor graph associated with the probabil- 
ity density ([!]). 



provided, so that the log-likelihood function for the model can be written as 



M 



(3) 



log L(s\g) = T - log Z(g) + 5 ^ 



M =i 



and depends on the data through the empirical averages 0^ = A X/t=i m Maximizing 
such function in order to find the coupling set g* best describing the data s in absence of 
prior leads to the set of M conditions 



(4) 



J P-I9 



where the quantities (0^) 9 = 9 log Z(g)/dg^ describe the averages of the functions 0^(s) 
under the model parametrized by g. Thus, in order to find the best model in order to 
describe data s one has to match empirical averages 0^ with ensemble averages {4>^) g *. 
An equivalent characterization of the problem of finding the optimal g* is provided by 
computing the Legendre transform 



(5) 



S((p) = max 
9 



M 



lo%Z(g) 



which can be identified as the Shannon entropy of the optimal model expressed as a func- 
tion of the empirical averages 0. Such quantity is the cumulant generating function for 
the maximum likelihood estimator of the couplings g* = <?*(0). From this perspective, the 
entropy S(4>) plays a role similar to the one of the free-energy F(g) = — log Z(g) in a statis- 
tical mechanics problem (i.e., estimating the ensemble averages {4>fj,) g given the couplings 



1 



IACOPO MASTROMATTEO 



g). In fact, F(g) satisfies the relations 

(6) (^)fl 



dF{g) 
d 2 F(g) 



dg^dg u 

being X/j,,u(g) the covariance matrix (cp^cp^g — {<fr[j,)g{4>u)g- Such matrix Xn,u is a central ob- 
ject in the field of information theory, where it is customarily referred as Fisher information 
[291 EQ]. For S$) it holds 

dS($) 



(8) 



(9) x^M 



d 2 s(4>) 



Just as the quantity Xn,u(g) can be used to express how fast in T the empirical averages 4> 
converge to the ensemble values {(ft), its matrix inverse expresses the rate of convergence of 
the inferred couplings g*(4>) to the asymptotic values g*(((fr}). In particular the covariance 
matrix for the estimator g*((fr) is given by 
(10) 

(g*M)-g;((<P)))(gt(4>)-gUm) -fim) (gUt)-gt((<P))) = %^ , 



It \ I t \ It 

where (-)t = Y[J=i p(s^). Hence the matrix x -1 (</0 expresses the stability of the in- 
ferred couplings with respect to the statistical error due to finite sampling. 

I will show in the following section how equations ^ and Q can be used to determine 
analytically g*{4>)- In particular, if it is possible to compactly express the entropy of a 
model as a function of the empirical averages the maximum likelihood estimator g*(4>) 
can be found easily. 

3. TWO FACTORIZABLE MODELS 

In this section I present two models for which the inverse problem can be exactly solved. 
The solution strategy can be straightforwardly extended to more general cases. 

In order to discuss this approach to the inverse problem, it is nevertheless necessary to 
introduce the notion of marginal pr, which given a set r C V is defined as 

(11) PT {sv) = Y, P^ s ) 

Si\i<£T 

and allows to define the cluster entropy Sr as 

(12) S r (pr) = - J>r(sr) logPr(sr) • 
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It is easy to prove that marginals and cluster entropies can be exactly expressed as functions 
of the averages of all the possible monomials contained in the cluster T. More precisely, 
given the set of 2' r ' monomials {</>r'}r'cr associated with the clusters V C r, it is easy to 
show that 

(!3) Pr(sr) = i £ 4>r'(sr>)(4>r>) g ■ 

r'cr 

This formula characterizes the marginals pr(sr) as local objects, as they are uniquely 
determined by the averages of potentials contained inside T. 



3.1. Tree-like graphs. Consider a model of the form 



(14) 



p(s) 



1 



Z(h, J) 



cxp 



where E = {(i, j) € V x V \ i < j} is a given set of edges. Suppose additionally that the 
graph associated with the edge set E is acyclic, so that the model (14) corresponds to an 



Ising model on a tree (more precisely, a forest). In that case, it can be shown by induction 
(see [301 02]) that; the total probability density p(s) can be factorized according to 



(15) 



p(s) = Yl P{ij}(8i,8j) Yl [P{i}(*0] 



1-1*1 



where di = {(j, k) £ E \ i £ (j, fc)}. Then, the entropy of the model is additive and can be 
written as 



(16) 



S(P)= ^2 S {i,j}( m i' m ^ C ^) + ^2( 1 -\ di \) S {i} 

1 V^T ft) 



mi 



where one has defined the magnetizations rrii = ^ Et=i s i anc ^ ^he correlations Cij = 

1 T (t) (t) i 1 

rEt=i*i s j • After the use of formula p3h, one finds by differentiation that the solution 



of the inverse Ising model on an acyclic graph is given by 

1 



(17) 



jEdi Si,Sj 

+ |(i-N)£ Si io g 



(1 + mjSj + rn-jSj + CijSiSj) 



(1 + mjSi) 



4 = !E s ^ io s 



-(1 + rriiSi + rn^Sj + c-ijSiSj) 
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and the stability of the solution is determined by 



(18) *^>'™ = 
The quantities 

(19) P{i}(sj) = ^(1 + rriiSi) 

(20) P{ij}( s i) s j) = j (1 + WjSj + + CijSiSj) 

appearing in the previous expression describe the empirical frequencies associated with the 
clusters {i} and {i,j}- 

Notice that the problem of finding any of the mean magnetizations (s^r^n or a correla- 
tions usually requires the use of iterative algorithms such as belief propagation, 
which converge in a number of steps linear with the number of vertices |V| [30J. For the 
inverse problem no iterative algorithm is required, as the solution can be found by simply 



evaluating formula (17). This is a consequence of the fact that - as observed in [271 EE] - 
inverse problems tend to be more local than their direct counterparts. Equivalently, due 
to sparsity of x -1 ; by shifting an empirical average (either m% or Jij) just neighboring 
couplings are changed, whereas for the direct problem the change of a coupling (say, hi), 
has an effect on a finite number of vertices (being the size of the perturbation roughly 
determined by the correlation length of the system) . These results are consistent with the 
ones of |20|, |3"T] , where the solution of the inverse Ising problem under the Bethe approxi- 



mation is derived. In particular, equation (17) can be recovered by specializing the results 



in \20\ [3"T] to the case of a tree, for which the Bethe approximation is exact. 

This solution requires indeed the knowledge of the full topology specified by E, so that this 
approach cannot be used to determine the graph structure from data. Rather, it allows to 
exactly recover the coupling strength given an acyclic set of edges E. 

3.2. One-dimensional models. Within this approach it is also possible to find the exact 
solution of the inverse problem for a one-dimensional chain with interactions of arbitrary 
order. Even if a heuristic solution for this model was first proposed in [32j, its detailed 
derivation is provided in the following. The model is defined by a probability distribution 
p{s) on a set of binary spins s E { — 1, 1} specified by a family of potentials acting on the 
first R spins <p(s\, . . . , sr) = (</>i(si, • • ■ , sr), . . . , 4>m{si, • • • , sr)) subject to the periodic 
boundary conditions Sj = Si + ^. By using translation operators T = {T n }^Q \ defined 
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through their action on the 4> 

(21) T n M (>i, ...S R ) = ^(Sl+np, ■ • • , SR+np) , 

one can define a one-dimensional chain as the probability distribution 



(22) 



p(s) 



1 / M N/p-1 \ 



Suppose that the family (f> = (fa, . . . , 4>m) contains all and only the 2^(1 — 2~ p ) monomials 
describing the unit cell (i.e., any potential of range smaller than R along the chain can be 
uniquely obtained by translating one of the 4>p by mean of a T n operator). Then one can 
exploit the factorization 



(23) 



p(s) = n ^l 



where one has defined the sets T n = {np+l, . . . , np+R} and 7 n = {(n+l)p+l, . . . , np+R}. 
This leads to the additivity formula for the entropy 

(24) S(p) = Y J [Sr n (Pr n )-S ln ( Pln )} 

n 

proved in appendix |Aj Translational symmetry and formula ( 13 ) lead finally to 



(25) 



E 

S 70 



E 



i M 

15 Yl E c ^ r ' <^<M s rO 



2 R 



m (1 r'cr M=i 
1 M 



2 R- P 



log 



1 M 

r'cr m=i 
1 A/ 

7'£70 M=l 



where the 4>r(sr) are monomials used to expand the marginals as in (13), while Cu,r = 1 if 



it exists an n such that T, 



nYp, 



br and c„ r = otherwise. The couplings are given by 



(26) 



E 

E 



,K E c /^,r"^r'(sr') 



2* 



r'cr 







log 



1 M 



2^ 



r'cr n v=\ 



TO 



Cjtt,7' 0y (sy) 



7'C7q 



log 



1 M 

j— 5^<V,y0f 0y(ay' 



2 R-p 
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while the inverse susceptibilities can be written as 

Sr'cr c mT' ^r(sr') Sr'cr c ^,r' <Ar'( s r' 



-l P 

X 



W 1 £- Pr.W) 



(27) - E 



[p S 7 'C 7o C M,7' (MV) P S 7 'C 70 C ^,7' 7 '( s 7') 



The structure of this solution is analogous to the one of ( |17[ ), in the sense that once that 
the S(p) is written as a sum of cluster entropies Sr n (pr n ) and S^ n (p ln ), the expression for 
the couplings is obtained by summing their separate contributions to the estimator g^(T(p). 
The main difference in this case is that one can exploit translational symmetry in order to 
equate the cluster contributions due to different unit cells. Thus, the number of sampled 
unit cells is y^TN/p, rather that \/T. This rules the convergence of the couplings to their 
asymptotic values: if the system is large translational symmetry allows to obtain reliable 
estimations of the couplings even with a small number of observations. 

4. Applications 

In order to validate the results shown above, some examples involving synthetic datasets 
are discussed. In the case of a tree-like model, a set of T = 10 6 configurations of a 
systems of size N = 50 has been simulated. The configurations have been obtained by 
MonteCarlo sampling from a distribution p(s) specified by a random tree-like graph E, 
while the couplings J and h have been uniformly drawn from the interval [0,1]. The 
variances of the inferred couplings for specific instances of the problem have been shown in 
figure[2j where they are compared against the expected scaling 1/T. The plot indicates that 



formula (17) correctly predicts the value of the inferred couplings and of their fluctuations. 
I have also considered the problem in which every coupling is multiplied by an inverse 
temperature parameter (3, in order to control the overall intensity of the fluctuations. The 
presence of strong noise corresponds to the regime of large /3, while the one in which 
noise is suppressed is associated with the small /3 region. In this example the couplings J 
and h have been uniformly drawn in [0, 1] and multiplied by a (3 parameter in the range 
[0.01,10]. Subsequently, T = 10 5 samples have been extracted via MonteCarlo. In figure 
[3] the variance of the inferred couplings is plotted against the inverse temperature f3. One 
can notice from the plots that, as the inverse problem depends just upon the products /3J 
and f3h, it is not possible to distinguish between an overall strength of the fluctuations 
a the temperature parameter f3. Hence, while the lowest inference error for the products 
f3J and fih is obtained for hi = Jy = 0, the lowest error for the values J and h (given a 
fixed, known value for (3) is obtained by finding a compromise between maximum signal 
(favoring high couplings) and minimum noise (favoring high temperature, or equivalently 
low j3). One also finds that the quality of the reconstruction degrades progressively in the 
large /3 regime by increasing the inverse temperature, until eventually some of the two-state 
marginals are not sufficiently sampled, in which case the inferred couplings diverge. The 
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Figure 2 . Variance of the inferred couplings h (red line) and J (blue line) 
against the number of samples T for a pairwise tree. The dashed green line 
plotted for reference indicates the expected scaling 1/T. 
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Figure 3. Variance of the inferred couplings h* (red line) and J* (blue 
line) against the inverse temperature /3 for a pairwise tree, obtained by 
using T = 10 5 MonteCarlo samples. Both the variances of h* and J* (left 
panel) and the ones of the products h*/3 and J*f3 (right panel) are plotted. 
The dashed lines indicate the estimations of the error obtained by using 



expression (18). 



same type of study has been carried on for the model (22) 
results obtained for a system of size N 
and periodicity parameter p = 2 



As an example, I show the 
50, characterized by interactions of range R = 4 
A number of configurations T = 10 6 have been simulated 



via MonteCarlo for a system with random interactions (g^ uniformly extracted in [0, 0.01]). 
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The results are shown in figure |4j where the variance of the inferred couplings is plotted 
against the expected scaling p/NT. Also for the one dimensional case I analyzed the effect 
of an inverse temperature parameter (3 £ [0.01, 1] modulating an interaction strength 
randomly drawn from the interval [0,1]. Figure [5] shows the errors in the reconstructed 
couplings as a function of the f3 parameter. Notice that the error on f3g* is minimum 




10 6 



Figure 4. Variance of the inferred coupling vector g* (blue line) plotted 
against the number of sampled unit cells NT/p, obtained by MonteCarlo 
sampling of a model describing a complete one-dimensional periodic chain 
of size N = 50, range R = 4 and periodicity p = 2. The green dashed line 
shows the reference scaling p/NT. 

for a flat distribution (/3 = 0), while the one on g* is minimized by an intermediate value 
of (5. The quality of the reconstruction is determined by how uniformly is the i?-spin 
marginal sampled. If local fluctuations are not sampled well-enough, the error on the 
inferred couplings is large as predicted by equation (27). As observed above, it is not 



necessary to probe global fluctuations of the system in order to accurately reconstruct the 
couplings. 

5. General structure of exact solutions 



What has been shown for models (14) and (22) can be easily extended to more general 
scenarios. In particular any factorization property for the probability distribution p(s) 



analogous to (15) and (23) breaks the entropy in a set of cluster contributions of the form 
(28) 



S(p) 



Teg 



c r Sr(pi 



where Q is a suitably chosen collection of clusters associated with a set of {crjreS co ~ 
efficients. Equation (13) implies that cluster entropies Sr{pr) are local in the empirical 
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Figure 5. Variance of the inferred couplings g* (blue line) against the 
inverse temperature (3 for a one dimensional periodic chain. T = 10 5 con- 
figurations have been sampled via MonteCarlo so to construct the empirical 
averages <f>. Both the the results for the inferred couplings g* (left panel), 
and the one for the products /3g* (right panel) are shown. The dashed line 



indicates the estimation of the error obtained by using equation (27) 



observables associated with T. Hence, if the model is large enough (i.e., if <p includes all 
the monomials located within the clusters T £ Q), the cluster entropies can be expressed 
as functions of the observables 4> conjugated to the couplings g. In this case the following 
facts hold: 

• Due to linearity of the derivative any inferred coupling g^(4>) receives contributions 
just by clusters containing its conjugated potential 0^(s). 

• The matrix x _1 is sparse, so that by shifting any of the empirical observables ^ 
the only clusters which are affected are the ones including variables contained inside 
the monomial <?^(s). 

The expressions (18) and (27) for the stability matrix x have then a clear interpretation 
in term of cluster contributions: if the clusters necessary to describe p{s) are uniformly 
sampled (i.e., the empirical frequencies pr are approximately flat), the error on the inferred 
couplings is small. Conversely, as soon as one ore more configurations sr are not sampled, 
the asymptotic error diverges. This is consistent with the observation that if a pr(sr) 
assigns zero mass to some configurations, infinite couplings are required to describe the 
marginals (in that case either some regularization can be used to enforce finiteness of the 
solution, or infinite solutions may be accepted as describing a hard constraint on the set 
of accessible configurations). As a further comment, notice that a reliable estimation of 
the couplings g doesn't require sampling global fluctuations in the configuration space: it 
is simply necessary to sample fluctuations which are local with respect to clusters in Q in 
order to accurately estimate the couplings g. 



These ideas can be useful even in the case in which the entropy doesn't strictly break 
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in the sum of small cluster terms. In particular it has been observed in \27\ [28] that for 
several models it is possible to account for a large fraction of the entropy S(p) by using 
just a few local terms (corresponding to the remark that in many cases the matrix x 1 is 
approximately sparse). In this case fast, although approximate, inference schemes could 
be obtained by using these techniques. 

Notice that the price to pay in order to analytically solve this inverse problem is the 
introduction of a potentially high number of parameters: a coefficient has to be fitted for 
each of the monomials included in any of the clusters r £ Q. In particular, the number of 
parameters to be considered is bound by M < \G\ X 2 maxree l r L Thus, even if a number of 
clusters polynomial in N is usually sufficient to properly describe the entropy of a system, 
their size enters exponentially in the number of fitted parameters. 



6. Conclusions 

The task of determining a probability distribution whose associated momenta match 
a given set of values is a known hard problem, with relevant implications in the field 
of high-dimensional inference. In this work I have shown that its exact solution can be 
found whenever the probability density for the system has a factorized structure. The 
sparsity of the Fisher information matrix x (0) an d the additivity of the entropy can 
be shown to rely upon this independence property. These result is completely general (in 
particular, it is independent of the interaction structure of the system). In order to provide 
illustrative examples, these ideas have been used to solve two specific inverse problems: 
the inverse Ising model on a tree and the one-dimensional periodic chain with arbitrary 
order interactions. The interactions of both systems can be accurately reconstructed if 
the clusters describing the system are uniformly sampled, while the inferred couplings 
diverge as soon as some relevant micro states are unobserved. These findings have been 
confirmed by numerical simulations, which show an excellent agreement with the analytical 
predictions. Techniques for the approximate solution of generic inverse problems could in 
principle be derived by developing the methods presented in this work. 



Appendix A. Factorization property for the one-dimensional chain 



I show in the following that for a one-dimensional periodic chain defined as in (22), the 
factorization property 

N/p-l 



(29) p(s)= ][ 



PT„(3 in ) 

„_„ Pm(» 7 ") 
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holds, where T n = {np+ 1, . . . ,np + R} and 7„ = {(n + l)p + 1, . . . , np + R}. To obtain this 
result, one can consider the auxiliary two-dimensional model defined by the log-probability 

N/p-l 

\ogp X (s,t) = - log Z(g)+ E Ev^( s l+np>---> s £+n P ) 

n=0 p£<f> 
N/p-l np+R 

(so) + a e E [(<? - o 2 + (<? - *r f1 ) 2 ] , 

n=0 i=(n+l)p+l 

in which the configuration space contains the original degrees of freedom are sf S { — 1,1} 
(with n = 0, . . . , N/p — 1 and i = 1 + np, ...,!? + rap) and the auxiliary ones t" S {—1, 1} 
(with n = 0, . . . , N/ p — 1 and i = 1 + (n + l)p, • • • , R + rap). The relation between the 
original model and the auxiliary one is sketched in figure [6} In particular the coupling A 
controls the strength of the bond in the auxiliary dimension (labeled by n) , so that the limit 
A — > oo describes the original chain with the obvious identification sf — > Sj and —> Sj. 
By defining the row variables s n = {s"}^Z^!^f and t n = {tf}l=f+(l?+i) p i ^ ne log-probability 



R 

■* ► 

n — 1 ##### s 1 

V" • • • K 
n = 2 m • • • • r 

• • • f 

n = 3 

• • • f 

Figure 6. Two dimensional auxiliary model p x (s) associated with the orig 
inal distribution p(s) describing a one-dimensional periodic chain. 



for the two dimensional model can be written as 

N/p-l 

(31) logp A 0M) = -logZ A ( 5 )- E 

n=0 

Hence the distribution over the degrees of freedom s n and t n defines a tree, because only 
successive row of variables interaclrl The marginals associated with the clusters T n and 7 n 



%nU n ) +^ti,n(s n if) + 7-(-n,n+l(t n ' ^ ) 



Periodic boundary conditions enforce the presence of a single loop of length N, so that the model is 
not exactly a tree. Nevertheless, for large enough and for g sufficiently distant from critical points of the 
model, if any, the presence of such loop can be neglected. 
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can be used in order to express the probability p (s,t) as 
^ p x, s f s = nn_Pr n u 7 n^ ?1 '^ ri )^7nur n+1 U ri ^ ?1+1 ' 



UnPrM n )Pi^ n ) 

where for the two-dimensional model T n and j n are analogously defined. By taking the 
A — > oo limit, the identification 

(33) p£ nU (s n ,t n ) ► pr n (snp+i,..-,s np+ R) 

A— ^oo 

(34) PrA^) ~ ► Pr n (Snp+l,---,S n p+R) 

A— >oo 



(35) P-y n (f) > Py n {s {n+1)p+1 ,...,s np+R ) . 

A— >oo 

allows to recover the factorization property ( |23| ). 
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